#ifndef Config_H
#define Config_H


#include <vector>
#include <string>
#include <map>
#include <queue>

#include "rapidjson/document.h"
#include "rapidjson/writer.h"
#include "rapidjson/stringbuffer.h"
#include "mfem.hpp"
#include "exprtk/exprtk.hpp"



//=================global variables: declare ===============================
extern const double smallValue;
extern const double verySmallValue;
extern const double largeValue;
extern const double veryLargeValue;

//length***Value is set in Mesh.cpp
extern double lengthSmallValue;
extern double lengthVerySmallValue;
extern double lengthLargeValue;
extern double lengthVeryLargeValue;
//=================global variables===============================


#include "Log/Log.h"
//#include "Time/SimulationTime.h"

class SimulationTime;
class IO;
class ParallelMPI;
class SpeciesCollection;
class RecordCollection;
class ParallelMesh;
class ElectroMagneticField;

namespace EasyMesh
{
    class Mesh;
    class ParallelMesh;
}

class Material
{
public:

    std::string name;

    std::string kind;
    
    double resistance;

    double heatCapacity;

    double conductivity;

    double relativePermittivity;

    double thermalConductivity;

};

class MeshParameter
{
public:
    std::string meshName;

    double meshScale;

    std::vector<std::string> mesh_periodic_edge0;
    std::vector<std::string> mesh_periodic_edge1;
    std::map < std::string, std::string > periodicFaceGroupMap;
    std::map < std::string, std::string > periodicFaceGroupMapReverse;

    std::map<std::string, std::string> physicalGroupToMaterialMap;
};

class FieldBoundaryConditionParameter
{
public:
    std::string field_name;

    //boundary name, like: left, right, inlet, outlet...
    std::string fbc_name;

    //boundary type: dirichlet, neumann, robin, float, open, none, externalCircuit
    std::string fbc_type;

    //for dirichlet, float, and neumann
    double fbc_value;

    std::string fbcValueString;

    double fbc_robin_k;

    double fbc_robin_gd;

    double fbc_robin_gn;


    //===========================================================================begin==================================================================================
    //accumlatedCharge(金属边界累计沉积电荷)，absorbedCurrent（边界吸收电流）等变量按理说应该在SpeciesCollection类中定义，但是这些变量在电势边界条件计算过程中经常用到，
    //因此为了统一管理，将这些变量都定义在场边界条件中。
    //
    //悬浮电势边界方法1： 
    //for float boundary condition: aCoef + bCoef*FloatPhi = Q(accumlatedCharge)
    //ref: PEGASUS Theoretical Manual page12
    //accumlatedCharge: 也作为悬浮电极上的初始电荷
    double accumlatedCharge;
    double aCoef;
    double bCoef;

    //悬浮电势边界方法2: 电容矩阵法
    //dVi = Pij * dQj, 其中dQj就是absorbedCharge，当前时间步增加的电荷量
    double absorbedCharge;
    double absorbedChargeInAccelerateStep;


    //absorbed current, can be used as power control
    //存储不同边界面的电流，用于放电功率控制，以及电流历史记录
    double absorbedCurrent;
    double absorbedCurrentAverage;
    std::queue<double> absorbedCurrentQueue;
    double fixedPower;
    int averageStep;

    //isDeposit = true表示是介质表面，电荷沉积加速的时候需要累计上面的电荷，只要有粒子边界条件是deposit就需要把对应边界的isDeposit设置为true
    //对应的场边界条件通常设置为none
    bool isDeposit;

    double capacityCii;

    //金属电极所受电场力，由麦克斯韦应力张量方法计算
    double electricForce[3];


    //FieldBoundaryConditionParameter中的累积变量（accumlatedCharge）需要在solve之前由setBoundaryCurrent设定，必须在Species::move, ParticleSource Psi,之后设定，因为这三个都会改变absorbedCharge
    //另外，在Species::move之前将absorbedCharge和absorbedCurrent重置为0，不能EM::solve之后，因为record还要输出边界电流absorbedCurrent
    //===========================================================================end==================================================================================

    //=======================================================externalCircuit begin====================================================================================
    //externalCircuitType: "RLC", 
    std::string externalCircuitType;

    //eca, ecb: externalCircuit coef a and coef b, 对应开发者文档式(3-16-11)或者式(3-16-21)
    double eca;
    double ecb;

    //externalCircuitType = "RLC"
    double rlcR;
    double rlcL;
    double rlcC;

    //rlcQCt: t时间步电容C上的电荷；rlcQCt1：t-dt时间步电容C上的电荷，也用于存储初始时刻电容C上的电荷；rlcQCt2：t-2dt时间步电容C上的电荷。
    double rlcQCt;
    double rlcQCt1;
    double rlcQCt2;

    double rlcKn;
    double rlcAlpha0;
    double rlcAlpha1;
    double rlcAlpha2;


    //externalCircuitType = "SWC", "SWC" is shunt two capacitor, 并联两个电容器，参考开发者文档[3.16.4 并联两个电容外电路]
    double swcC1;
    double swcC2;

    //=======================================================externalCircuit end======================================================================================

    FieldBoundaryConditionParameter() 
    {
        absorbedCharge = 0.0;
        accumlatedCharge = 0.0;
        absorbedCurrent = 0.0;
        absorbedCurrentAverage = 0.0;
        averageStep = 10;
        fixedPower = 0.0;
        isDeposit = false;
        absorbedChargeInAccelerateStep = 0.0;
    };
    FieldBoundaryConditionParameter(FieldBoundaryConditionParameter& fieldBoundaryConditionParameter) 
    {
        absorbedCharge = 0.0;
        accumlatedCharge = 0.0;
        absorbedCurrent = 0.0;
        absorbedCurrentAverage = 0.0;
        averageStep = 10;
        fixedPower = 0.0;
        isDeposit = false;
        absorbedChargeInAccelerateStep = 0.0;
    };


    void compile()
    {
        //init  for exprtk
        symbol_table.add_variable("t", t);

        expression.register_symbol_table(symbol_table);

        parser.compile(fbcValueString, expression);
    }

    //fbc_value 和calValue()有时候不用的，在外电路边界中，calValue()是时变电源的电压值，而fbc_value是电极上的电压值，这两个通常不一样
    double calValue(double tInp)
    {
        t = tInp;

        double result = 0.0;
        result = expression.value();
        return result;
    };

    exprtk::symbol_table<double> symbol_table;
    exprtk::expression<double> expression;
    exprtk::parser<double> parser;

    double t;
};


class ResistanceConnectionParameter
{
public:
    std::string boundaryName1;
    
    std::string boundaryName2;
    
    double resistance;
};

class ElectroMagneticFieldParameter
{
public:

    std::string externalFieldString;

    //是否自洽求解电磁场，默认：ture
    bool isSolveEM;

    //define poisson equation solver: hypre, petsc, SuperLU_serial, SuperLU_mpi, 默认：hypre
    std::string solverType;

    double solverRelativeTolerance;
    int solverNumSweeps;
    int solverMaxLevels;
    int solverRelaxType;
    int solverRelaxOrder;

    std::string externalMagneticFieldFromPegasusFileName;

    std::string externalMagneticFieldString;


    std::vector<FieldBoundaryConditionParameter*> fbc_param;


    //====================== 电器连接：电压连接，电容矩阵法计算悬浮电势所用到的参数 ===========================================================
    
    //electricalConnectionVector: 电器连接组，每组包含一个或多个金属边界名称（字符串），输入文件中只包含大于等于2个边界的组，求解器中应将单个悬浮边界加入到组中
    //只考虑"float"悬浮边界
    std::vector< std::vector<std::string> > voltageConnectionVector;

    //floatBoundaryMethod，悬浮电势边界求解方法：iterationMethod, capacitanceMatrixMethod
    std::string floatBoundaryMethod;

    //====================== end ==============================================================================================================


    //====================== 电器连接：电阻连接 =================================================================================================

    //electricalConnectionVector: 电器连接组，每组包含一个或多个金属边界名称（字符串），输入文件中只包含大于等于2个边界的组，求解器中应将单个悬浮边界加入到组中
    //只考虑"float"悬浮边界
    std::vector<ResistanceConnectionParameter> resistanceConnectionParameterList;

    //====================== end ===============================================================================================================

};


class ElectroMagneticSourceParameter
{
public:

    //sourceType: totalScatterdField, 
    std::string sourceType;

    std::string boundaryName;

    std::vector<double> waveInitialPosition;

    std::vector<double> totalFieldRegionCenterPosition;

    std::vector<double> totalFieldRegionSize;

    double electricFieldAmplitude;

    //运动方向 与 z轴 的夹角（默认 90°）
    double angleTheta;

    //运动方向在x-y平面投影 与 x轴 的夹角（默认 0°）
    double anglePhi;

    //极化角， 即 在 θ̂ -ϕ̂  平面上， 电场方向 与 θ̂  的夹角（默认 180°）
    double polarizationAngleAlpha;

    //pulseType: "sin", "gaussian"
    std::string pulseType;

    double gaussianPulseTimeWidth;

    double sinPulseFrequency;

    double plasmaInnerRadius;

    double plasmaThickness;

};

class ElectroMagneticBoundaryConditionParameter
{
public:


    std::string boundaryName;

    //boundaryConditionType: PEC，extrapolating, none, 
    std::string boundaryConditionType;

};


class RcsParameter
{
public:

    //rcsType: "monostatic", "bistatic"
    std::string rcsType;

    double farFieldRadius;

    double anglePhi;

    double angleTheta;

    double gaussianFrequencyMax;

    int gaussianFrequencyNumber;
};


class ElectroMagneticFieldDGTDparameter
{
public:

    //有限元的阶数
    int femOrder;

    //默认是TM波； 否则就是TE波
    bool isTM;

    //磁导率，默认国际单位制
    double magneticPermeability;

    //介电常数，默认国际单位制
    double dielectricConstant;

    //完美匹配层相关参数
    double upmlThickness;
    int upmlM;
    int upmlN;
    bool isReal;

    //总场边界发射的是否是sin波，还是高斯波Exp
    //bool isSin;
    //bool isExp;

    bool isCalRCS;
    //每隔extr_steps步，保存一次总场散射场上的数据
    int extr_steps;

    std::string electronDensityString;

    std::string electronDensityFileName;


    std::vector<ElectroMagneticSourceParameter> electroMagneticSourceParameterList;

    std::vector<ElectroMagneticBoundaryConditionParameter> electroMagneticBoundaryConditionParameterList;

    RcsParameter rcsParameter;
};





class SimulationTimeParameter
{
public:

    //time step for the simulation
    double dt;

    //number of total timesteps to perform in the simulation
    int timestep_simulation;

    int timestep_io;

    int timestep_average;

    //timestep for all kind collision
    //timestep for DSMC and coulomb should be several times of this timestep
    int timestep_collision;

    //number of timestep of performing coulomb collision
    int timestep_coulomb;

    //number of timestep of performing DSMC collision
    int timestep_DSMC;

    //number of timestep of performing DSMC collision
    int timestep_restore;

    bool isReload;

    int reloadIOstep;



    //=============== 卫星表面电荷沉积加速参数：begin===========================================
    bool isChargeDepositAccelerate;
    int timestepSteady;
    int timestepAccelerate;
    double maxAccelerateFactorK;
    //=============== 卫星表面电荷沉积加速参数：end=============================================

    


};




class ParticleBoundaryConditionParameter
{
public:

    //edge name, like: left, right, inlet, outlet...
    std::string edge_name;

    //particle boundary type: absorb, reflect, periodic
    std::string pbc_type;

    std::string periodic_edge_name;
};


class SpeciesParameter
{
public:

    std::string species_name;

    int species_number;
    
    //species_type: particle, fluid
    std::string species_type;
    
    //move_type: neutral, e(only electric field), em(electric field and magnetic field)
    std::string move_type;

    std::string init_position_type;

    std::string init_velocity_type;

    //coefficient  of particle max number
    double c_part_max;

    double mass;

    double charge;

    double number_density;

    std::string numberDensityString;

    //n_part_per_cell_for_weight is used to calculate particle weight
    //Now weight is the same for all the particles!!!
    //if initialize paritcles before PIC loop, the n_part_per_cell should be equal to n_part_per_cell_for_weight
    int n_part_per_cell;
    //int n_part_per_cell_for_weight;
    int maxMacroParticleNumberPerCell;
    //如果粒子数大于maxMaxMacroParticleNumberPerCell，即使粒子权重大于combineMaxWeight仍继续合并，为了处理电离倍增问题
    int maxMaxMacroParticleNumberPerCell;
    double combineMaxWeight;
    double combineMaxPhysicalParticleNumber;
    int combineMaxMacroParticleNumberPerCell;
    std::string combineCriterion;
    std::string combineAlgorithm;

    //all particles of one species has the same weight
    //weight = nb_density / n_part_per_cell; if n_part_per_cell = 0, weight = nb_density / 100
    double weight;


    //temperature, isotropic
    double temperature;

    //temperature, anisotropic
    double temperature_x;
    double temperature_y;
    double temperature_z;

    std::vector<double> mean_velocity;

    //control timesteps, use larger timesteps for some heavy particles
    int timestep_frozen;

    //Boundary conditions for particules
    std::string bc_part_type_west;
    std::string bc_part_type_east;
    std::string bc_part_type_south;
    std::string bc_part_type_north;
    std::string bc_part_type_bottom;
    std::string bc_part_type_up;


    //=============== parameters for DSMC ===================================
    double diameter;
    //reference temperature
    double ref_temperature;
    double visc_temp_index;
    double vss_scat_inv;

    //==============Parameters for physical sputtering and backscattering============
    //atomic number: H 1; D 1; C 6; Ar 18;
    int atomic_number;

    //atomic mass: unit is amu, H 1.008; D 2.016; T 3.024; Be 9.012182; Ar 39.95
    double atomic_mass;

    //target surface_binding_energy for physical sputtering, W: 11.75(old value 8.68) C: 7.41
    double surface_binding_energy;

    //target density, unit: g/cm^3 （Ai = 10^-10 m）, for physical sputtering, C: 2.248 g/cm^3, W: 19.35 g/cm^3
    double density_solid;

    //============Parameters for backscattering as the target
    //number of constituent elements in the target.
    int ne2;

    //array for atomic numbers of the constituents.
    std::vector<int> an2_vector;

    //array for relative numbers of the constituents.
    std::vector<int> nw2_vector;

    //=============particle boundary condition==============
    std::vector<ParticleBoundaryConditionParameter> pbc_param;
    std::map<std::string, std::string> particleBoundaryTypeMap;

    //是否输出粒子信息
    bool isWrite;

};


class ParticleSourceParameter
{
public:

    std::string species_name;

    std::string species_name_dependent;

    int species_number;

    int species_number_dependent;

    //particle source type: load or emit
    std::string ps_type;

    std::vector<double> timeDuration;

    int timestep_update;

    bool is_dependent;

    //false: add particle source only at first timestep; true: add particle source every timestep
    bool is_every_time;

    //==================================================Load particle source=======================================
    std::string load_type;

    double number_density;

    //load particle number is fixed, load time step number is calculated
    double load_particle_number;

    //load density per second
    double dn;

    double dn_scale_down;
    
    // dq = dn (m-3) * temperature (eV)
    double dq;

    double temperature;

    std::vector<double> mean_velocity;

    std::vector<double> load_position_start;
    std::vector<double> load_position_end;

    //box[x0, x1, y0, y1]
    std::vector<double> load_box;

    std::string profileFilterString;

    //isFIxNTV = true: fix plasma density, temperature, and velocity in the source region
    bool isFixNTV;

    //==================================================Load particle source: end =======================================
       
    //emitType: flux, rate, current density, current, thermalElectronEmit, driftMaxwellian, photoElectricEmission, fluxVaryingWithTime
    std::string emitType;

    //emitCoordinateSystem: local(发射定向速度沿表面法线,emit_velocity指定), global(发射定向速度基于全局三维坐标，emitVelocityVector指定)
    std::string emitCoordinateSystem;

    double emit_particle_number;

    double emit_flux;
    std::string emitFluxString;

    std::string emit_position;

    double emit_offset;
    
    std::string edge_name;
    
    double emit_velocity;

    std::vector<double> emitVelocityVector;

    bool isDepositOppositeCharge;

    //emitVelocityType: directionMaxwellian, driftMaxwellian, 目前默认是directionMaxwellian，对用户不提供设置选项。
    std::string emitVelocityType;

    //parameters for thermalElectronEmit: Richardson-Dushman model
    //ref: 《空间电推进》-于达仁 P238
    double workFunction;        //unit: eV
    double surfaceTemperature;         //unit: K, the temperature of emission object surface
    double A;                   //Richardson coefficient, constant, 120A/(〖cm〗^2∙K^2)



    //parameters for photoElectricEmission
    //ref: 开发者文档3.10 光照产生光电子处理方法
    double emitCurrentDensity;      //发射电流密度   
    double scatteringAngle;         //散射角
    std::vector<double> lightDirection;     //光照方向



};


class CollisionParameter
{
public:

    std::string reaction_name;

    std::string reaction_type;
   
    std::vector< std::vector<double> > cross_section;

    std::vector<std::string> species_name_list;
    std::vector<std::string> species_left_list;
    std::vector<std::string> species_right_list;

    std::vector<int> species_number_list;

    int reaction_id;

    std::string collisionAlgorithm;

    int timestep_collision;

    double time_zoom_factor;

    std::string cs_file_name;

    //========================= parameters for field ionization: begin ============================================
    double ionizationEnergy;
    //========================= parameters for field ionization: end ============================================
};

class PSIParameter
{
public:

    std::string edge_name;

    std::string psi_type;

    std::vector<std::string> species_name_list;

    std::vector<int> species_number_list;

    //position offset of injected or sputtered particles
    double emit_offset;

    //if is_self_consistent is false, only calculate the psiRate (like sputtering rate or reflection rate)
    //not create new particles
    bool is_self_consistent;

    //for 1d: left or right
    std::string psi_position;

    //for recycling
    double recycling_factor;

    double recycling_temperature;

    double temperature;

    //for sputerring

    //for secondaryEmit

    //secondary emit coefficient
    std::string sec;

    double randomSecLimit;

    int emitPlusNumber;

    //velocityType: maxwell(half maxwell, out of surface), singleEnergy(enery is single, angle is random), beam(one direction normal to surface), cos
    std::string velocityType;

    std::string psiAlgorithm;

    bool isDepositOppositeCharge;


};



class DiagnosticParameter
{
public:

    //vdf: velocity distribution funtion
    //edf: energy distribution function
    //vedf: some parameters of vdf and edf are same

    //division number of velocity and energy
    int vedf_n_velocity;
    int vedf_n_position;
    int vedf_n_species;

    std::vector<double>         vedf_position;
    std::vector<std::string>    vedf_species_name_list;
    std::vector<int>            vedf_species_number_list;

    std::vector<double>         edf_max_energy;

};


class RecordParameter
{
public:

    //recordType: particleMacroNumber, particlePhysicalNumber, fieldAtMonitorPoint
    std::string recordType;

    std::string recordName;

    std::vector<double> monitorPointPosition;


    //boundary name for cal RecordThruster
    std::string faceGroupName;
};


class LicenseParameter
{
public:

    //licenseVersion: singleMachine; network
    std::string licenseVersion;

    //std::string serverIp;
    //std::string serverPort;
};



//---------------------------------------------------------------------------------------------------------------------
//Config class: holds all the properties of the simulation that are read from the json input file
//---------------------------------------------------------------------------------------------------------------------
class Config 
{

 public:

    Config(ParallelMPI* pmpiIn);


    void read_species();

    void read_particle_source();

    void read_collision();

    void read_psi();

    void read_diagnostic();

    void read_simulation_time();

    void readMeshParameter();

    void readElectroMagneticFieldParameter();

    void readElectroMagneticFieldDGTDparameter();

    void readMaterialDatabase();

    void readRecord();

    void readLicense();

    void compute();

    void add_simulation_time(SimulationTime* sim_time_);

    void addParallelMPI(ParallelMPI* pmpiInp);

    int get_species_number(const std::string& species_name);

    bool get_string(std::string& cpp_value, std::string value_name);
    bool get_double(double& cpp_value, std::string value_name);
    bool get_int(int& cpp_value, std::string value_name);
    bool get_bool(bool& cpp_value, std::string value_name);

    bool get_string_vector(std::vector<std::string>& cpp_value, std::string value_name);
    bool get_double_vector(std::vector<double>& cpp_value, std::string value_name);
    bool get_int_vector(std::vector<int>& cpp_value, std::string value_name);

    bool get_class_string(const rapidjson::Value& json_value, std::string& cpp_value, std::string value_name);
    bool get_class_double(const rapidjson::Value& json_value, double& cpp_value, std::string value_name);
    bool get_class_int(const rapidjson::Value& json_value, int& cpp_value, std::string value_name);
    bool get_class_bool(const rapidjson::Value& json_value, bool& cpp_value, std::string value_name);


    bool get_class_string_vector(const rapidjson::Value& json_value, std::vector<std::string>& cpp_value, std::string value_name);
    bool get_class_double_vector(const rapidjson::Value& json_value, std::vector<double>& cpp_value, std::string value_name);
    bool get_class_int_vector(const rapidjson::Value& json_value, std::vector<int>& cpp_value, std::string value_name);

    static bool greater_sort_pbc(ParticleBoundaryConditionParameter& e1, ParticleBoundaryConditionParameter& e2)
    {
        return (e1.edge_name < e2.edge_name);
    }

    //=========== json related variables ========================
    const char* file_name = "input/easypic3d.json";

    rapidjson::Document dom;


    //========== Config memeber variables for simulation =====
 
    int dimension;

    //electromagneticFieldType: "electrostatic", "electromagnetic", "none"
    std::string electromagneticFieldType;

    //only used for 2 dimension
    bool isCylindricalCoordinateSystem;

    //有限元求解电磁场的阶数
    int order;


    //bool isReload;
    
    bool isBarycentricCoordinateParticleTracking;
    
    bool isWriteAllProcess;

    bool isWriteVtk;

    std::string dataBaseDir;

    int debugLevel;

    //用于控制是否启用介质中电荷移动的功能，主要用于卫星表面充电的介质漏电现象,只有3D支持此功能
    bool isChargeInDielectricMove;

    double zDirectionLength;


    std::vector<SpeciesParameter> species_param;

    std::vector<ParticleSourceParameter> particle_source_param;

    std::vector<CollisionParameter> collision_param;

    std::vector<PSIParameter> psi_param;

    std::vector<DiagnosticParameter> diagnostic_param;

    SimulationTimeParameter sim_time_param;

    LicenseParameter license_param;

    SimulationTime* sim_time;

    MeshParameter mesh_param;

    ElectroMagneticFieldParameter em_param;

    ElectroMagneticFieldDGTDparameter emDGTD_param;

    std::map<std::string, Material> materialMap;

    std::vector<RecordParameter> record_param;


    //some physics constants
    const double const_c;
    const double const_e;
    const double const_emass;
    const double const_ephi0;
    const double const_ephi0_inv;
    const double const_pi;
    const double const_boltz;     //Boltzmann constant
    const double const_h;         //Planck constant
    const double const_unitAtomMass;


    //global class points
    IO* io;
    EasyMesh::ParallelMesh* mesh;
    ParallelMPI* pmpi;
    SpeciesCollection* speciesCollection;
    RecordCollection* recordCollection;
    ElectroMagneticField* emField;

    //database
    std::map<std::string, Material*> dataBaseMaterialMap;

};

#endif
